{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Continuous Ranked Probability Score (CRPS) in probabilistic forecasting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In point estimate forecasting, the model outputs a single value that ideally represents the most likely value of the time series at future steps. In this scenario, the quality of the predictions can be assessed by comparing the predicted value with the true value of the series. Examples of metrics used for this purpose are the Mean Absolute Error (MAE) and the Root Mean Squared Error (RMSE).\n",
    "\n",
    "In probabilistic forecasting, however, the model does not produce a single value, but rather a representation of the entire distribution of possible predicted values. In practice, this is often represented by a sample of the underlying distribution (e.g. 50 possible predicted values) or by specific quantiles that capture most of the information in the distribution.\n",
    "\n",
    "One of the main applications of probabilistic forecasting is the estimation of prediction intervals - ranges within which the actual value is expected to fall with a certain probability. In this case, the model should aim to achieve the desired coverage (e.g. 80%) while minimising the width of the prediction interval.\n",
    "\n",
    "The Continuous Ranked Probability Score (CRPS) is a generalisation of the Mean Absolute Error (MAE) tailored to probabilistic forecasting. Unlike the MAE, which compares point predictions to observations, the CRPS evaluates the accuracy of an entire predicted probability distribution against the observed value. It does this by comparing the empirical cumulative distribution function (CDF) of the predicted values with the step-function CDF of the true value. \n",
    "\n",
    "Two key components of the CRPS are the empirical CDF of the predicted values, 𝐹(𝑦), and the CDF of the observed value, 𝐻(𝑦). The CRPS is then calculated as the integral of the squared difference between these two functions over the entire real line:\n",
    "\n",
    "+ Empirical CDF of the forecast, $F(y)$: This is constructed from the ensemble of predicted values. Each predicted value contributes a \"step\" in the cumulative distribution. The predicted values are therefore treated  as a sample of the underlying probability distribution.\n",
    "\n",
    "+ CDF of the observed Value, $H(y)$: This is a step function that transitions from 0 to 1 at the true observed value. It represents the probability that the observed value falls below a given threshold.\n",
    "\n",
    "The CRPS measures the area between the two CDFs, $F(y)$ and $H(y)$, across all possible values of $y$. Mathematically, it is expressed as:\n",
    "\n",
    "\n",
    "$$\\text{CRPS}(F, H) = \\int_{-\\infty}^{\\infty} \\big(F(y) - H(y)\\big)^2 \\, dy$$\n",
    "\n",
    "\n",
    "This integral quantifies the squared difference between the forecasted and observed distributions. \n",
    "\n",
    "The CRPS can be computed for a single observation or for a set of observations. In the latter case, the CRPS is averaged over all observations to provide a summary measure of the model's performance.\n",
    "\n",
    "CRPS is widely used in probabilistic forecasting because it provides a unified framework for evaluating both the sharpness (narrowness) and calibration (accuracy) of predictive distributions. By doing so, it ensures that models are not only accurate in their point predictions but also appropriately represent uncertainty. Smaller values of CRPS indicate a better match between the forecast and the observed outcome."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style=\"text-align: center\">\n",
    "<img src=\"../img/crps.gif\" style=\"width: 600px;\">\n",
    "<br>\n",
    "<font size=\"2.5\"> <i>Example of CRPS calculation</i> </font>\n",
    "</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CRPS and Skforecast"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Skforecast provides different output options for probabilistic forecasting, two of which are:\n",
    "\n",
    "- **`predict_bootstrapping`**: Returns multiple predicted values for each forecasted step. Each value is a variation of the forecast generated through bootstrapping. For a given step \\(i\\), \\(n\\) predictions are estimated.\n",
    "\n",
    "- **`predict_quantile`**: Returns the estimated values for multiple quantiles. Internally, the forecaster uses `predict_bootstrapping` and then calculates the desired quantiles.\n",
    "\n",
    "For both outputs, the CRPS (Continuous Ranked Probability Score) can be calculated to evaluate the forecasting performance of the model.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"admonition note\" name=\"html-admonition\" style=\"background: rgba(0,191,191,.1); padding-top: 0px; padding-bottom: 6px; border-radius: 8px; border-left: 8px solid #00bfa5; border-color: #00bfa5; padding-left: 10px; padding-right: 10px;\">\n",
    "\n",
    "<p class=\"title\">\n",
    "    <i style=\"font-size: 18px; color:#00bfa5;\"></i>\n",
    "    <b style=\"color: #00bfa5;\">&#128161 Tip</b>\n",
    "</p>\n",
    "\n",
    "<p>For more examples on how to use probabilistic forecasting, check out the following articles:</p>\n",
    "<ul>\n",
    "    <li>\n",
    "        <a href=\"https://cienciadedatos.net/documentos/py42-probabilistic-forecasting\" target=\"_blank\">\n",
    "            Probabilistic forecasting with machine learning\n",
    "        </a>\n",
    "    </li>\n",
    "    <li>\n",
    "        <a href=\"https://cienciadedatos.net/documentos/py60-probabilistic-forecasting-prediction-intervals-multi-step-forecasting\" target=\"_blank\">\n",
    "            Probabilistic forecasting: prediction intervals for multi-step time series forecasting\n",
    "        </a>\n",
    "    </li>\n",
    "    <li>\n",
    "        <a href=\"../faq/probabilistic-forecasting-crps-score.html\" target=\"_blank\">\n",
    "            Continuous Ranked Probability Score (CRPS) in probabilistic forecasting\n",
    "        </a>\n",
    "    </li>\n",
    "</ul>\n",
    "\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CRPS from a sample of predictions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Continuous Ranked Probability Score (CRPS) is calculated by comparing the empirical cumulative distribution function (ECDF) of the forecasted values to the step function CDF of the true value. When the available information consists of the true value (`y_true`) and a sample of predictions (`y_pred`), the CRPS can be calculated by following these steps:\n",
    "\n",
    "1. Generate the Empirical Cumulative Distribution Function (ECDF) of the predictions:\n",
    "   - Sort the predictions.\n",
    "   - Use each sorted prediction as a step in the ECDF.\n",
    "\n",
    "2. Generate the Cumulative Distribution Function (CDF) of the true value:\n",
    "   - Since there is only a single true value, this is represented as a step function that jumps from 0 to 1 at the observed value (`y_true`).\n",
    "\n",
    "3. Calculate the CRPS by integrating the area between both curves:\n",
    "   - Create a grid of values to evaluate the ECDF. This grid is the combination of the predictions and the true value.\n",
    "   - Compute the squared differences between the forecasted ECDF and the true CDF, and then summing the areas between the two curves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/joaquin/miniconda3/envs/skforecast_19_py13/lib/python3.13/site-packages/pymc_extras/model/marginal/graph_analysis.py:10: FutureWarning: `pytensor.graph.basic.io_toposort` was moved to `pytensor.graph.traversal.io_toposort`. Calling it from the old location will fail in a future release.\n",
      "  from pytensor.graph.basic import io_toposort\n"
     ]
    }
   ],
   "source": [
    "# Libraries\n",
    "# ======================================================================================\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from skforecast.metrics import crps_from_predictions\n",
    "from skforecast.metrics import crps_from_quantiles\n",
    "from scipy.interpolate import interp1d\n",
    "import properscoring as ps\n",
    "from CRPS import CRPS\n",
    "from pymc_marketing.metrics import crps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "np.float64(2.021484883451464)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Simulate data: true value and and array of 100 predicted values for the same true value\n",
    "# ======================================================================================\n",
    "y_true = 500\n",
    "y_pred = np.random.normal(500, 10, 100)\n",
    "crps_from_predictions(y_true, y_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result of <code>skforecast.metrics.crps_from_predictions</code> is compared with other implemented functions in the `properscoring`, `CRPS` and `pymc_marketing` libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "properscoring : 2.0214848834514623\n",
      "CRPS          : 2.021484883451464\n",
      "pymc-marketing: 2.021484883451463\n"
     ]
    }
   ],
   "source": [
    "# properscoring, CRPS, pymc-mar\n",
    "# ==============================================================================\n",
    "print(f\"properscoring : {ps.crps_ensemble(y_true, y_pred)}\")\n",
    "print(f\"CRPS          : {CRPS(y_pred, y_true).compute()[0]}\")\n",
    "print(f\"pymc-marketing: {crps(y_true, y_pred.reshape(-1, 1))}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When forecasting multiple steps, the CRPS can be calculated for each step and then averaged to provide a summary measure of the model's performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>y_true</th>\n",
       "      <th>y_pred</th>\n",
       "      <th>crps_from_predictions</th>\n",
       "      <th>properscoring</th>\n",
       "      <th>CRPS</th>\n",
       "      <th>pymc_marqueting</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>90.108786</td>\n",
       "      <td>[8.640637789445714, -1.3080015845984816, 12.14...</td>\n",
       "      <td>82.100538</td>\n",
       "      <td>82.100538</td>\n",
       "      <td>82.100538</td>\n",
       "      <td>82.100538</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>96.322133</td>\n",
       "      <td>[3.1678558225107745, 6.737363230274925, 5.6735...</td>\n",
       "      <td>88.410864</td>\n",
       "      <td>88.410864</td>\n",
       "      <td>88.410864</td>\n",
       "      <td>88.410864</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>112.879253</td>\n",
       "      <td>[6.709160916434245, 10.896201858093296, 0.9120...</td>\n",
       "      <td>105.460630</td>\n",
       "      <td>105.460630</td>\n",
       "      <td>105.460630</td>\n",
       "      <td>105.460630</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>101.939744</td>\n",
       "      <td>[14.521434285699028, 1.1295876122380442, 15.13...</td>\n",
       "      <td>94.259885</td>\n",
       "      <td>94.259885</td>\n",
       "      <td>94.259885</td>\n",
       "      <td>94.259885</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>109.202309</td>\n",
       "      <td>[13.80532539228533, 5.482203757147254, 6.46324...</td>\n",
       "      <td>101.908526</td>\n",
       "      <td>101.908526</td>\n",
       "      <td>101.908526</td>\n",
       "      <td>101.908526</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       y_true                                             y_pred  \\\n",
       "0   90.108786  [8.640637789445714, -1.3080015845984816, 12.14...   \n",
       "1   96.322133  [3.1678558225107745, 6.737363230274925, 5.6735...   \n",
       "2  112.879253  [6.709160916434245, 10.896201858093296, 0.9120...   \n",
       "3  101.939744  [14.521434285699028, 1.1295876122380442, 15.13...   \n",
       "4  109.202309  [13.80532539228533, 5.482203757147254, 6.46324...   \n",
       "\n",
       "   crps_from_predictions  properscoring        CRPS  pymc_marqueting  \n",
       "0              82.100538      82.100538   82.100538        82.100538  \n",
       "1              88.410864      88.410864   88.410864        88.410864  \n",
       "2             105.460630     105.460630  105.460630       105.460630  \n",
       "3              94.259885      94.259885   94.259885        94.259885  \n",
       "4             101.908526     101.908526  101.908526       101.908526  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# CRPS for multiple steps\n",
    "# ==============================================================================\n",
    "rng = np.random.default_rng(123)\n",
    "n_steps = 40\n",
    "n_bootstraps = 100\n",
    "predictions = pd.DataFrame({\n",
    "    'y_true': rng.normal(100, 10, n_steps),\n",
    "    'y_pred': [rng.normal(5, 5, n_bootstraps) for _ in range(n_steps)]\n",
    "})\n",
    "\n",
    "predictions['crps_from_predictions'] = predictions.apply(lambda x: crps_from_predictions(x['y_true'], x['y_pred']), axis=1)\n",
    "predictions['properscoring'] = predictions.apply(lambda x: ps.crps_ensemble(x['y_true'], x['y_pred']), axis=1)\n",
    "predictions['CRPS'] = predictions.apply(lambda x: CRPS(x['y_pred'], x['y_true']).compute()[0], axis=1)\n",
    "predictions['pymc_marqueting'] = predictions.apply(lambda x: crps(x['y_true'], x['y_pred'].reshape(-1, 1)), axis=1)\n",
    "display(predictions.head())\n",
    "\n",
    "assert np.allclose(predictions['properscoring'], predictions['CRPS'])\n",
    "assert np.allclose(predictions['properscoring'], predictions['pymc_marqueting'])\n",
    "assert np.allclose(predictions['crps_from_predictions'], predictions['properscoring'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "crps_from_predictions    93.974753\n",
       "properscoring            93.974753\n",
       "CRPS                     93.974753\n",
       "pymc_marqueting          93.974753\n",
       "dtype: float64"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Average CRPS\n",
    "# ==============================================================================\n",
    "mean_crps = predictions[['crps_from_predictions', 'properscoring', 'CRPS', 'pymc_marqueting']].mean()\n",
    "mean_crps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CRPS from quantiles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A quantile is a value that divides a data set into intervals, with a specific percentage of the data lying below it. Essentially, it is a point on the cumulative distribution function (CDF) that represents a threshold at which a given proportion of the data is less than or equal to that value.\n",
    "\n",
    "For example, the 40th percentile (or 0.4 quantile) is the value below which 40% of the data points fall. To find it, you would examine the CDF, which shows the cumulative proportion of the data as you move along the values of the data set. The 0.4 quantile corresponds to the point where the CDF reaches 0.4 on the vertical axis, indicating that 40% of the data lies at or below this value.\n",
    "\n",
    "This relationship between quantiles and the CDF means that, given several quantile values, it is possible to reconstruct the CDF. This is essential for calculating the **Continuous Ranked Probability Score (CRPS)**, which measures the accuracy of probabilistic forecasts by comparing how well the predicted distribution matches the true value.\n",
    "\n",
    "Given a set of quantiles, their associated probabilities, and the true value, the CRPS can be calculated as follows:\n",
    "\n",
    "1. Construct the Empirical Cumulative Distribution Function (ECDF) using the quantiles and their corresponding probabilities.\n",
    "   \n",
    "2. Generate the CDF for the true value: Since the true value is a single point, its CDF is represented as a step function that jumps from 0 to 1 at the observed value.\n",
    "   \n",
    "3. Calculate the CRPS as the squared diference between the two curves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "np.float64(1.7339183102042313)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# CRPS score from quantiles\n",
    "# ==============================================================================\n",
    "y_true = 3.0\n",
    "\n",
    "quantile_levels = np.array([\n",
    "    0.00, 0.025, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55,\n",
    "    0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, 0.975, 1.00\n",
    "])\n",
    "pred_quantiles = np.array([\n",
    "    0.1, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5,\n",
    "    8.0, 8.5, 9.0, 9.5, 10.0, 10.5, 11.0, 11.5\n",
    "])\n",
    "\n",
    "crps_from_quantiles(y_true, pred_quantiles, quantile_levels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, results are compared versus the `properscoring` package. In this case, a warapper function is used to calculate the CRPS score from the predicted quantiles using `crps_quadrature`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "np.float64(1.7342500001706027)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def crps_from_quantiles_properscoring(y_true, predicted_quantiles, quantile_levels):\n",
    "    \"\"\"\n",
    "    Calculate the Continuous Ranked Probability Score (CRPS) for a given true value\n",
    "    and predicted quantiles using the function crps_quadrature from the properscoring\n",
    "    library.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    y_true : float\n",
    "        The true value of the random variable.\n",
    "    predicted_quantiles : np.array\n",
    "        The predicted quantile values.\n",
    "    quantile_levels : np.array\n",
    "        The quantile levels corresponding to the predicted quantiles.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    float\n",
    "        The CRPS score.\n",
    "    \"\"\"\n",
    "    if len(predicted_quantiles) != len(quantile_levels):\n",
    "        raise ValueError(\n",
    "            \"The number of predicted quantiles and quantile levels must be equal.\"\n",
    "        )\n",
    "    \n",
    "    # Ensure predicted_quantiles are sorted\n",
    "    sort_idx = np.argsort(predicted_quantiles)\n",
    "    predicted_quantiles = predicted_quantiles[sort_idx]\n",
    "    quantile_levels = quantile_levels[sort_idx]\n",
    "\n",
    "    def empirical_cdf(x):\n",
    "        # Interpolate between quantile levels and quantile values\n",
    "        cdf_func = interp1d(\n",
    "            predicted_quantiles,\n",
    "            quantile_levels,\n",
    "            bounds_error=False,\n",
    "            fill_value=(0.0, 1.0),\n",
    "        )\n",
    "        return cdf_func(x)\n",
    "\n",
    "    # Integration bounds\n",
    "    xmin = np.min(predicted_quantiles) * 0.9\n",
    "    xmax = np.max(predicted_quantiles) * 1.1\n",
    "\n",
    "    # Compute CRPS\n",
    "    crps = ps.crps_quadrature(np.array([y_true]), empirical_cdf, xmin, xmax)\n",
    "\n",
    "    return crps[0]\n",
    "\n",
    "\n",
    "crps_from_quantiles_properscoring(y_true, pred_quantiles, quantile_levels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Results are similar but not identical. This may be due to differences in the implementation of the CRPS calculation or the numerical methods used to approximate the integral. The skforecast team is working on validating the implementation of the CRPS function in the library."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "skforecast_19_py13",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.9"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
